home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qk.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-12-09  |  3.2 KB  |  103 lines
  1. /* integration/qk.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <float.h>
  22. #include <math.h>
  23. #include <gsl/gsl_integration.h>
  24. #include "err.c"
  25.  
  26. void
  27. gsl_integration_qk (const int n, 
  28.                     const double xgk[], const double wg[], const double wgk[],
  29.                     double fv1[], double fv2[],
  30.                     const gsl_function * f, double a, double b,
  31.                     double *result, double *abserr,
  32.                     double *resabs, double *resasc)
  33. {
  34.  
  35.   const double center = 0.5 * (a + b);
  36.   const double half_length = 0.5 * (b - a);
  37.   const double abs_half_length = fabs (half_length);
  38.   const double f_center = GSL_FN_EVAL (f, center);
  39.  
  40.   double result_gauss = 0;
  41.   double result_kronrod = f_center * wgk[n - 1];
  42.  
  43.   double result_abs = fabs (result_kronrod);
  44.   double result_asc = 0;
  45.   double mean = 0, err = 0;
  46.  
  47.   int j;
  48.  
  49.   if (n % 2 == 0)
  50.     {
  51.       result_gauss = f_center * wg[n / 2 - 1];
  52.     }
  53.  
  54.   for (j = 0; j < (n - 1) / 2; j++)
  55.     {
  56.       const int jtw = j * 2 + 1;    /* j=1,2,3 jtw=2,4,6 */
  57.       const double abscissa = half_length * xgk[jtw];
  58.       const double fval1 = GSL_FN_EVAL (f, center - abscissa);
  59.       const double fval2 = GSL_FN_EVAL (f, center + abscissa);
  60.       const double fsum = fval1 + fval2;
  61.       fv1[jtw] = fval1;
  62.       fv2[jtw] = fval2;
  63.       result_gauss += wg[j] * fsum;
  64.       result_kronrod += wgk[jtw] * fsum;
  65.       result_abs += wgk[jtw] * (fabs (fval1) + fabs (fval2));
  66.     }
  67.  
  68.   for (j = 0; j < n / 2; j++)
  69.     {
  70.       int jtwm1 = j * 2;
  71.       const double abscissa = half_length * xgk[jtwm1];
  72.       const double fval1 = GSL_FN_EVAL (f, center - abscissa);
  73.       const double fval2 = GSL_FN_EVAL (f, center + abscissa);
  74.       fv1[jtwm1] = fval1;
  75.       fv2[jtwm1] = fval2;
  76.       result_kronrod += wgk[jtwm1] * (fval1 + fval2);
  77.       result_abs += wgk[jtwm1] * (fabs (fval1) + fabs (fval2));
  78.     };
  79.  
  80.   mean = result_kronrod * 0.5;
  81.  
  82.   result_asc = wgk[n - 1] * fabs (f_center - mean);
  83.  
  84.   for (j = 0; j < n - 1; j++)
  85.     {
  86.       result_asc += wgk[j] * (fabs (fv1[j] - mean) + fabs (fv2[j] - mean));
  87.     }
  88.  
  89.   /* scale by the width of the integration region */
  90.  
  91.   err = (result_kronrod - result_gauss) * half_length;
  92.  
  93.   result_kronrod *= half_length;
  94.   result_abs *= abs_half_length;
  95.   result_asc *= abs_half_length;
  96.  
  97.   *result = result_kronrod;
  98.   *resabs = result_abs;
  99.   *resasc = result_asc;
  100.   *abserr = rescale_error (err, result_abs, result_asc);
  101.  
  102. }
  103.